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ABSTRACT 


The class of one-dimensional stretching functions used in finite- 
difference calculations is studied. For solutions containing a highly 
localized region of rapid variation, simple criteria for a stretching 
function are derived using a truncation error analysis. These criteria 
are used to investigate two types of stretching functions. One is an 
interior stretching function, for which the location and slope of an 
interior clustering region are specified. The simplest such function 
satisfying the criteria is found to be one based on the inverse 
hyperbolic sine. It was first employed by Thomas et al. The other 
type of function is a two-sided stretching function, for which the 
arbitrary slopes at the two ends of the one-dimensional interval are 
specified. The simplestsuch general function is found to be one 
based on the inverse tangent. The special case where the slopes were 
both equal and greater than one was first employed by Roberts. The 
general two-sided function has many applications in the construction of 
finite-difference grids. An example of such an application is found in one 
of the references. 



I. INTRODUCTION 


Finite-difference calculations of fluid flow problems are best carried 
out using an equispaced grid in a rectangular (or cubic) computational 
domain, with the flow variables and components of the position vector 
as dependent variables, and boundary conditions applied at the edges 
(or faces) of the domain. In order to minimize the number of grid 
points required for a given accuracy, one seeks boundary-fitted coordinate 
transformations that cluster points in regions where the dependent 
variables undergo rapid variation. These regions may be due to body 
geometry (very large curvatures or corners), compressibil ity (entropy 
layers, shock waves and contact discontinuities), and viscosity (boundary 
layers and shear layers). A complex flow may thus contain a variety 
of such regions of various length scales, and often of unknown location. 

An ideal grid would adjust with each time or iteration step to maintain 
optimum clustering. Such adaptive grid methods, which involve the 
solution of auxiliary equations, have been developed for one-dimensional 
problems (refs. 1-3). Their extension to complex multi -dimensional 
flows is a difficult problem, particularly when the regions requiring 
clustering do not have simple topological properties required by a 
finite-difference grid. 

There are many practical problems in which the locations and length 
scales of regions of rapid variation can be estimated a priori 
(e.g., known geometry, attached boundary layers, simple shock wave 
configurations) . In these cases the clustering can be incorporated 
in automatic grid generators which solve an elliptic boundary-value 
problem (refs. 4-6). The distribution of grid points on the boundaries 
is then normally prescribed algebraically, using one-dimensional 
stretching functions. (Here, stretching function refers to any 
transformation involving stretching or clustering). It is also possible 
to employ stretching functions to obtain a clustered grid from an 
unclustered grid by applying clustering to one coordinate family only. 

For some simple geometries, one can construct entire clustered grids 
purely algebraically, using only one-dimensional stretching functions. 
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The simplest class of one-dimensional stretching functions is that 
involving two parameters. In interior stretching functions, the 
parameters are the location and slope of a single clustering region. 

In two-sided stretching functions, the slopes at the two ends of the 
one-dimensional interval are specified. The anti -symmetric two-sided 
stretching function (with the same slope at each end) is of special 
interest, since the portion from the midpoint to either end defines 
a one-sided stretching function (with zero curvature at the other 
end). A one-sided stretching function can also be obtained as a special 
case of the interior stretching function, by locating the clustering 
region at one end. Since the clustered end has zero curvature in this 
case, these two one-sided stretching functions are of a different nature. 

An interior stretching function based on the inverse hyperbolic sine was 
employed by Thomas, Vinokur, Bastianon, and Conti (ref. 7) in a 
numerical solution of inviscid supersonic flow over a blunt delta wing 
with elliptical cross section. The function was used to cluster points 
on the body at the vertices of very eccentric ellipses. The one-sided 
version clustered points in the flow near the body surface to resolve 
the entropy layer due to the bow shock. No derivation of the stretching 
function was given, and the clustering parameter appearing in it was not 
related to the length scales of regions of rapid variation in physical 
space. 

An anti -symmetric two-sided stretching function of a logarithmic type 
was employed by Roberts (ref. 8) to study boundary layer flows. The 
heuristic derivation of the function avoided consideration of the 
truncation errors associated with finite-difference approximations. 

While this function has been used successfully in many flow calculations, 
there is a need for a general two-sided function which allows arbitrary 
stretching or clustering to be specified independently at each end. An 
application would be problems in which the appropriate length scales 
requiring clustering are significantly different at the two ends. Another 
application is the distribution of grid points on a curve which is 
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defined piecewise, where continuity of grid spacing is desired at the 
ends of the piecewise segments. A third application is the use of such 
a function as a blending or interpolating function to construct two 
and three-dimensional grids using one-dimensional stretching functions 
and shearing transformation. The construction of a single, well- 
ordered grid for wing body flows by Vinokur (ref. 9) is an example of 
these applications. 

The present work has two objectives. One is to obtain simple, rational 
criteria for one-dimensional stretching functions, by considering the 
truncation errors inherent in f inite-difference approximations. The 
functions introduced by Thomas et al and Roberts will be found to be 
the simplest ones satisfying these criteria. The other objective is to 
derive a simple form of the general two-sided stretching function. 
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II. TRUNCATION ERROR ANALYSIS FOR ONE-DIMENSIONAL STRETCHING FUNCTIONS 


An exact analysis of the truncation error inherent in a finite-difference 
calculation would require knowledge of the equation being solved and the 
finite difference approximation that is used. Here we are concerned with 
the special situation where the solution contains a highly localized 
region of rapid variation with respect to some coordinate, and we seek 
approximate criteria for a stretching function that will be independent 
of the equation or difference algorithm. The quantities that are 
approximated are in general non-linear functions of the unknowns and 
their spatial derivatives. The error analysis will be performed in 
terms of the fractional truncation errors for the spatial derivatives. 

Let r(t) be the equation describing a ^-coordinate curve, where t is 
any parameter that varies smoothly with arc length. If A and B denote 
the ends of the curve, we introduce the normalized variables 

t = (t - t A )/(t B - t A ) 
and 5 = (f - ? A )/(c B - ? A ), 


(la) 

(lb) 


ranging from 0 to 1 . For simplicity, all partial derivatives with 
respect to 5 or t will be written as total derivatives. Let 4>( t) be 
any function of the unknowns. Outside of the region where d<j>/dt = 0, 
we can define a natural length scale of the variation of <p with respect 
to t as 



di> /aj> 

dt 2 


( 2 ) 


Since the components of r also enter as dependent variables in the 
calculation of metrics and Jacobians, we similarly define the natural 
length scale of the variation of r with respect to t as 


L 


rt 


dr 

1 / 

d r 

dt 

/ 

dt 2 


(3) 
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Note that if t is proportional to arc length, then L t is precisely 
the radius of curvature, normalized by the length of the curve, 


Assume first that t is used as the normalized computational variable 
(i.e. £ = t), with At as the uniform grid spacing. Let 5<j>/6t denote the 
finite difference approximation to d^/dt . Using equation (2), one can 
write any first order accurate approximation as 

ft-* -jjf-n ♦ou- 1 +tl ta. (4) 


Similarly, the approximation to dr/dt can be written as 


6r - [1 + 0 (L _1 


fit 


rt 


At)] 


(5) 


If L"^ t or L’^ rt became very large in some localized region, then a 
prohibitively small At would be required to obtain a desired fractional 
truncation error. Outside of the localized region, the excessive number 
of grid points would be wasted. The obvious remedy is to seek a new 


-1 


-1 


computational variable t, for which L ' and L" 1 remain of 0(1), even 
-1 -1 $?> ^ 


though L ^ or L rt could be locally very large. 


With the aid of the identities 


d<j> _ dj> dt 
d£ ” dt d? 

and 

d 2 4> _ d 2 <j> /dt_>2 d^t 

dt 2 ' dt 2 ^ dt dt 2 - 


and definition (2), one can easily show that 


-1 < , -1 

' « § 


+ L 


-1 


t5 , 


( 6 ) 


(7) 


( 8 ) 
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Similarly, using definition (3), one obtains the inequality 


L 1 s l - 1 dt_ + l“1 
rc rt dc tc- 

One criterion for the stretching function £(t) is therefore 

L '' te ■ 0 o. 


( 9 ) 


( 10 ) 


In addition, we require that 


L 'V§ = ° 0 )> 


and 


-1 


dt 


L _ ‘ rt dc = 0(1) 


(ID 


( 12 ) 


Consider the case where L~^ t is very large in some localized region. It 
follows from equation (11) that 


i =0 V 


(13) 


in that region. Since L~\. remains large over an interval whose 
thickness is of 0(L,.), we require that dt/d£ remains small over that 

'pt _1 

interval. Furthermore, since L ^ = 0(1) outside of the localized 
region, it follows from equation (11) that dt/d? = 0(1), i.e. that dt/d£ 
does not become large anywhere. But these two additional requirements 
are satisfied if condition (10) is valid. Noting that 


d /dt - _ d t //dt v 

dt W d? 2 / d ? ’ 

it follows upon integration over a finite interval At that 


(14) 


a /dt , 

4 ( w) 


'• < L \>max 4t 


(15) 
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Applying equation (15) to the localized region, for which At = 0(L t ), 
we find, using equation (10), that 


. /dt. 


= 0 ^t l 


(16) 


Thus equation (13) is satisfied over the entire localized region. Letting 
At = 1 in equation (15), we see that dt/d£ = 0(1) is satisfied over the 
complete range of t. In the event that L ^ is larger than L ^ in the 
localized region, then condition (13) is replaced by 


dt _ 

dc 
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The above analysis is easily extended to higher order finite-difference 
approximations, as well as the treatment of higher derivatives. In order 
to consider fractional truncation errors due to a second-order accurate 
approximation to d<j>/dt (and also for regions where d 4>/ dt =0), it is 
appropriate to define a different length scale of the variation of <p with 
respect to t as 



We now require that L~l. remains of 0(1), even though L could be locally 

<p£ <pt 

very large. Using the identity 


d 3 ij> _ d 3 <() ,dt, 3 , ,d 2 <t> dt d 2 t , d£ d 3 t 

dt 2 d C 2 dt d£ 3 • 


(19) 


we obtain the inequality 




dts2 

dc ; 


3 d*L + [_"2 
J L <j>t dc tc L tc. 


( 20 ) 


In addition to satisfying equation (10), the stretching function £(t) 
must satisfy 


L 


-1 


U 


0 ( 1 ). 


( 21 ) 
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Also, if L is larger that L -1 ^ in the localized region, condition (13) 
must be replaced by 


dt 


0 <V. 


( 22 ) 


2 2 

A first order finite difference approximation to d <j)/dt requires the 
introduction of the length scale 


r _ d 2 $ /d 3 <. 

» t_ dt 2 / dt 3 . 

Combining equations (2), (18), and (23), we see that 


(23) 


I 



/ 



(24) 


Thus conditions (10), (21), and (13) or (22) are sufficient to guarantee 
that L'J remains of 0(1). The length scale L t is also the appropriate 
one to use at a point where d<)>/dt = o. At such a point, using equations 
(7) and (19), we obtain the relation 


C" 1 i L 

<K 


-1 


<pt 


dt 

d? 


3 L 


t€. 


(25) 


The criteria for ?(t) is again equation (10), with equation (13) replaced 
by 


dt 

dC 


0 <(>t^ 


(26) 


if the localized region of rapid variation occurs around the point dg>/dt = 0. 
Conditions (22) and (26) are to be replaced by analogous ones containing 
L rt and if these are the more significant length scales. 

In summary, one first defines length scales appropriate to the difference 
approximation and the location of the region of rapid variation. The 
criteria for the stretching function £(t) can be stated as: 
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1. All the inverse length scales of the variation of t with respect to £ 

must be at most of order one throughout the range of t. 

2. The slope dt/d£ must be of the order of the minimum length scale of 

the variation of <t> or r with respect to t in the localized region of 

rapid variation. 

These criteria will insure that most of the grid points will be concentrated 
in the localized region of rapid variation, with a sufficient number of 
points left in the remainder of the domain. The criteria will be used to 
investigate the two-parameter stretching functionsof the next two sections. 

III. A GENERAL TWO-SIDED STRETCHING FUNCTION 

In this section we derive a general two-sided stretching function 

E,{ t; Sq, s^), where £ and t are normalized variables defined by equation (1), 

and the parameters s Q and s-j are dimensionless slopes defined as 

s o ■ ft <°> (27a) 

and 

s r I"). (27*)' 

In order to be useful for contructing finite-difference grids, the 
function must be monotonic, and satisfy conditions (10) and (21) even if 
s Q or s.| becomes very large. For a general range of applications, it 
would be desirable for the function to be simple, invertible, and to vary 
continuously over the complete ranges of Sq and Sy 

An attractive candidate for such a stretching function is a scaled portion 
of a single, universal function w(z). For a given Sq and s^ , the stretching 
function will be obtained by properly scaling the portion of the universal 
function from corresponding points Zq and z-j . An additional requirement 
is that the unnormalized function £(t) be independent of the designation 
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of a particular end as A or B. One can easily show that this restricts 
tne universal function to be odd, i.e. 

w(-z) = -w(z). 


The simplest odd, monotonic, invertible functions are sinz and tanz. 

Their hyperbolic relatives are produced by letting z be complex. The 
inverse functions are formed by associating z with either t or £. One 
can determine whether either universal function is suitable as a basis 
for a stretching function by applying conditions (10) -\nd (21) for 
very large Sq or s y Actually, the most extreme test occurs for the 
simpler anti -symmetric casesg = s^, which corresponded to Zq = -z^, with 
z being either real or pure imaginary. 

An evaluation of the anti -symmetric two-sided stretching functions 
obtained from sinz and tanz is carried out for the case Sq = s^ > 1 
in Appendix A. Only tanz produces a stretching function satisfying 
conditions (10) and (21), with the inverse length scales being 
logarithmically of 0(1). The stretching function £(t) is a scaled 
portion of the inverse hyperbolic tangent. Expressing the hyperbolic 
tangent as a logarithm, we obtain exactly the function derived by 
Roberts (ref. 8). It turns out that is a piecewise linear function 

of t, a property that was used by Roberts to define his function. This 
suggests a related stretching function, for which L~^ is a piecewise 
linear function of £. The corresponding universal function is the 
error function erfz. The associated stretching function is also 
analyzed in Appendix A, and found to satisfy conditions (10) and (21). 
However, it is not invertible, and has larger maximum inverse length 
scales than the former stretching function. 

On the basis of the above considerations, the universal function w = tanz 
will be used to obtain a stretching function for arbitrary s Q and s^. 
Introducing the ranges 

Az = z 1 - z Q (28) 
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and 


Aw = tanz.| - tanzQ , 


(29) 


we can define the normalized variables 


S = (z - z 0 )/Az 


(30) 


and 


t = (tanz - tanz Q )/Aw. 

The slope of the stretching function is then given by 


d£ = 

dt 




Using the trigonometric identity 

sin (z, - z Q ) 
tanz l - tanz 0 = cos z ~ cos — 


1 


0 , 


we find for the parameters s Q and s-j the relations 

sinAz cos z Q 
s 0 ” Az cos z^ 


(31) 


(32) 


(33) 


(34) 


and 


sin Az cos z 


s, = 


1 


1 Az cos z. 


This suggests introducing the new parameters 


B =/^7 


(35) 


(36) 


and 

A = A 75 !". 


(37) 


11 


The parameters A and B can then be expressed in terms of Zq and z-j 
as 


and 


D _ sin Az 

D “ 7 — 

Az 


cos z. 


A = 


cos z 


( 38 ) 


(39) 


Using the cosine sum identity, we can also write equation (39) as 


A = cos Az + tanz-j sin Az 


(40a) 


and 


1/A = cos Az - tanZg sin Az- 


(40b) 


For a given value of B, Az is obtained by solving equation (38). Equation 
(40) can then be solved to obtain Zq and Aw for a given value of A. The 
stretching function obtained from equations (30) and (31) can then be 
written as 


t = 


tan (£Az + z Q ) - tanz Q 
Aw 


(41) 


While equation (41) is a formal expression for the general stretching 
function, it cannot be used for calculations in its present form. Depending 
on the value of B, Az and Aw are either real or pure imaginary. For certain 
ranges of A and B, z Q can become complex. U^ing the tangent sum identity and 
equation (39), we can eliminate z Q from equation (41) and obtain instead 


t = 


tan gAz 

A sin Az + (1 - A cos Az) tan £Az 


(42) 


This can be further simplified by noting that A = 1 corresponds to the 
anti -symmetric solution which was analyzed in Appendix A. Let us denote 
this solution as u(j;). Setting A = 1 in equation (42), and using the 
tangent sum identity, we obtain 
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(43) 


u = 1/2 + 


tan [AzUz 1/2) 1 
2 tan (Az/2) 


In terms of u, equation (42) then takes the simple form 


t = u 
Z " A + (1 - A)u , 

which can be readily inverted as 

.. = t 

U " (1/A) + (1 - 1 / A) t . 

Note that both u(t) and its inverse can be obtained as scaled portions 
of a rectangular hyperbola. For each function, the slopes at the two 
ends are A and 1/A, respectively. Finally, we observe that for 
calculational purposes, equations (44) and (45) are well behaved in the 
neighborhood of A = 1 . 


(44) 


(45) 


We thus have the remarkable result that one essentially needs to know 
only the anti -symmetric stretching function for the geometric mean of 
the slopes S Q and . The square root of the ratio of those slopes 
determines an additional simple transformation which produces the 
desired stretching function. Since both equations (43) and (44) are 
invertible, the resultant stretching function is also invertible. The 
key trigonometric property making this result possible is that the 
tangent of a sum is a rational function of the individual tangents. By 
contrast, the sine of a sum involves the individual sines and cosines, 
and is not expressible as a rational function of sines alone. An analysis 
of a stretching function based on w = sinz has been carried out, but is 
not presented here. The parameters corresponding to B and A turn out to 
be the arithmetic mean and difference of the two slopes. A separation 
into two functions corresponding to equations (43) and (44) is not 
possible, and one must use the direct form corresponding to equations 
(41) and (42). For B > 1, the inversion involves the solution of a 
quadratic equation, and the sign of A must be tested in order to choose the 
appropriate root. It is indeed serendipity that the tangent function 
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dictated by truncation error considerations is also the much simpler one 
for constructing a general two-sided stretching function. 


The calculation of the anti -symmetric function depends on the size of B. 
If B > 1, it follows from equation (38) that z is imaginary and we obtain 
the results (previously obtained in Appendix A) 


B = 


sinh Ay 
Ay 


(46) 


and 


u = 1/2 + 


tanh [Ay ( T - 1/2)] 
2 tanh (Ay/2) 


The inversion of equation (47) yields 

r - 1/2 + tanh" 1 [(2u - ptanh (Ay/2)'l 
^ ' 2Ay 


(47) 


(48) 


Note that the hyperbolic tangent and its inverse can be expressed in 
terms of exponentials and a logarithm, respectively. 


For B < 1, Az is real, and the corresponding results are 


and 


B = 


sin Ax 
Ax , 


u = 1/2 + 


5 = 1/2 + 


tan [Ax ( 5 - 1/2)] 

2 tan (Ax/2) , 

tan -1 [ ( 2u-l ) tan (Ax/2)] 
2A X 


(49) 

(50) 

(51) 


When B is very near one, both of the above formulations break down, 
since ax and Ay approach zero. The appropriate expressions are obtained 
by expanding equations (49) and (50) in powers of Ax. To first order 
in B-l, one obtains 


u = £ [1 + 2(B - 1) ( 5-.5)(l -5 )] (52) 

and 

5 - u [1 - 2(B - 1 ) (u - .5) (1 - u)]. (53) 
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By scaling half of the above functions, one obtains the one-sided 
stretching functions with 5 q given at t = 0 and zero curvature at t = 1 . 

The results are: 


s o > 1 

s 0 = 

sinh 2Ay 
2Ay 

(54) 

and 

t = 

tanh [Ay(5 - 1_H_ 
tanh Ay , 

, . tanh" 1 [ ( t - 1 )tanh Ay] 
1 + Ay 

(55) 


5 = 

(56) 

V 1 

s o : 

sin 2Ax 
2A X 

(57) 


t = 

, , tan [Ax(£ - 1 )] 
tanAx , 

(58) 


5 = 

. t tan" 1 [ ( t - 1 )t an ax] 

Ax 

(59) 

s o = 1 

t = 

5 [1 - .5(s q - 1)(1 - 0(2 - 0] 

(60) 

and 

5 = 

t [1 + .5(s q - 1)(1 - t) (2 - t)]. 

(61) 

The two-sided 

stretching functions for B>1 and one sided 

stretching 

function for 

V 1 

require the inversion of the function 



y = 

sinh x/x. 

(62) 


An approximate analytic representation of the inverse function x = f(y), 
valid over the range of y required by a stretching function, is derived 
in Appendix B. The results are as follows: 

For y< 2.7829681 

? -3 (63) 

x = /6y (1 - .15y + . 057321 429y - .024907295y 

+ . 0077424461 y 4 - .001 07941 23y 5 ) , 
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where 



y 

= y - 1. 

(64) 

For y> 2.7829681 


X = V 

+ (1 + 1/v) log (2v) - .02041793 + .24902722w 
+ 1.9496443W 2 - 2.6294547w 3 + 8. 5679591 lw 4 , 

(65) 

where 

v = log y 

(66) 

and 

w = 1/y - .028527431. 

(67) 


The maximum magnitude of the fractional error in y as defined by equation 
(B.2) is .000267732 for 1 <y< 69.64. The magnitude of the error reaches 
.00083 at y = 120.5. These errors are small enough so that the grids 
constructed by the resultant stretching functions will exhibit slope 
discontinuities that are negligible within the accuracy of any practical 
finite-difference approximation. 

The two-sided stretching function for B<1 and one-sided stretching function 
for Sq < 1 require the inversion of the function 

y = sin x/x. (68) 

An approximate analytic representation of the inverse function is derived 
in Appendix C. The results are as follows: 

For y <.26938972 (69) 

x = tt[1 - y + y 2 - (1 + Ti 2 /6)y 3 + 6.794732y 4 

- 1 3 . 205501 y 5 + 11.726095y 6 ]. 

For . 26938972 < y < 1 (70) 

x = >/6y (1 + .15y + . 057321 429y 2 + .048774238y 3 

- . 053337753y 4 + . 0758451 34y 5 ), 

where (71) 

y = i - y. 
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The maximum magnitude of the fractional error in y defined by equation 
(C.2) is .00019717, which is small enough for numerical applications. 

While it is appealing to view the general two-sided stretching function 
as a distortion of an anti -symmetric stretching function via equations 
(43) and (44), there are advantages in looking at the more basic forms 
of the solution given by equations (41) and (42). If efficiency in terms 
of operations count is important, then the optimum form for the case 
B > 1 is derived from equation (42), as either 

+ = tanh SAy . , 

L sinh Ay + (1 - A cosh Ay ) 

or t * e 254y -1 (72b) 

e 25Ay () . Ae -Ay, ; Ae ay 

The most efficient form for the case B<1 is equation (41) with z replaced 
by x throughout. 

It is also instructive to study the general solution (41), and see how it 
reduces to different real representations for various ranges of A and B. 
This is carried out in Appendix D. It is easy to demonstrate that for 
B < 1 , the solution is a scaled portion of a tangent, while for B = 1 
it is a scaled portion of a rectangular hyperbola. For B>1 (and the 
corresponding Ay given by equation (46), the representation depends on 
the value of A. As shown in Appendix 0, if exp (-Ay)<A<exp (Ay), 
the solution is a scaled portion of a hyperbolic tangent. Outside of 
that range the corresponding function is the hyperbolic cotangent, 
while on the boundary it is the exponential function. 

When A = 1, the stretching function is anti -symmetric, and the solution 
I curve contains an inflection point. As A departs from one, the inflection 

point moves towards one end, and eventually could disappear. It is shown 
in Appendix D that an inflection point will be present if 1/cosh Ay < A 
f < cosh Ay for B > 1 , and cos Ax < A < 1 /cos A x for B < 1 . If B < 2/n, the 

solution must always contain an inflection point. Figure D.l is a log-log 
1 plot of the B vs A plane, showing the inflection point regions, as well 
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as the boundaries for the various real representations of the solution. 
Another way to study equation (41) is to follow the path of the solution 
in the complex z-plane. This is also carried out in Appendix D, with the 
results plotted in Figure D.2. The plot shows clearly the transition from 
one representation to another as the parameters A and B are varied. 

A typical set of curves for the two-sided stretching function t ( 5 ) is 
shown in figure 1 for Sq = 100, and for values of s^ ranging from.l to 
100. The figure shows the smooth transition from solutions containing 
an inflection point, for large s^ , to those without an inflection point 
for small s^ . Note that the high concentration of points at the two ends 
for s-j = 100 still leaves a sufficient number of points to resolve the 
central region. 


IV. A GENERAL INTERIOR STRETCHING FUNCTION 


In this section we derive a general interior stretching function 
£(t; s^ t.j ) , where s i is the dimensionless slope at the inflection 
point t.j , i .e. 


s i 


= 

dt 


(V 


(73) 


and 


(t.) = 0 (0<t < 1). 

dt 1 


(74) 


We limit our consideration to s. >1, which is the only case of practical 
interest. We again look for a function which is a scaled portion of an 
odd universal function w(z). As in Section III, the simple functions 
sinz and tanz are considered first, and conditions (10) and (21) are 
examined for the anti -symmetric case (t^ = 1/2) when s^ is very large. 

The evaluation is carried out in Appendix E, and sinz is found to produce 
an appropriate function with inverse length scales that are logarithmically 
of 0(1). The stretching function £(t) is a scaled portion of the inverse 
hyperbolic sine. 
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The general interior stretching function for arbitrary t. is readily 
obtained from the universal function w = sinz, by letting z = iy. In 
terms of the range Ay = y^ - y^, and the implicity defined £..= £(t..), 
the final result can be written as 

sinh [Ay(C-Cj)] 
sinh Ay 

The inverse function is 

K =5j + ^ sinh' 1 [(t/t. - 1) sinh ^Ay]. (76) 

The relation between^ and t i (for a given Ay) is obtained from equation 
(75) by setting t = £ = 1. The result can be written as 




1=1- cosh Ay + sinh Ay coth ifiy. 


(77) 


Expressing the inverse hyperbolic cotangent in terms of a logarithm, we 
can write the inverse of equation (77) as 

fl + t i (e Ay - 1)' 


5 = J. 

i 


2Ay 


log 


L - t, 0 - 


(78) 


Equations (76) and (78) are precisely the ones given by Thomas et al 
(ref. 7) 


If s.j and t.j are the given parameters, the corresponding value of Ay 
must be calculated. This can be done by differentiating equation (75) 
and substituting into equation (73). Using equation (77) to eliminate 
C . , we can write the result as 


1 


(s i t i Ay)‘ 


cosh Ay - 1 + 1/t.j 


-12 


sinh Ay 


(79) 


This is an implicit equation for Ay involving two independent parameters. 
If the interior point is not too close to either end, and the slope s^ 
is sufficiently large, one can obtain a simplification. Assuming that 
exp (-2Ay)< <1, we can approximate equation (79) as 


19 


( 80 ) 


2s i 




(1 - t.): 


sinh (Ay/2) 
Ay/ 2 


Equation (80) is in the form of equation (62), and one can use the 
approximate analytic inversion derived in Appendix B. 


The special case of an anti -symmetric solution is obtained by setting 
t.j = = 1/2. The results (derived in Appendix E) are 

t ■ w [, , (BD 


and e, = 1/2 + ^ sinh 1 [(2t-l )sinh(Ay/2)], 

(82) 

_ _ sinh (Ay/2) . 

where b i n/2 Ay/2 

(83) 

A one-sided stretching function is obtained by setting t^ 
The results can be written as 

= 5i “ 0. 

t - sinh UAy) 
sinh Ay 

(84) 

and F - 1 sinh _1 (t sinh Ay) , 

^ " Ay 

(85) 

. . . _ sinh Ay 

where s i ' s o ' Ay 

(86) 


It is interesting to compare the one-sided stretching function derived 
from the hyperbolic tangent (equations (54) to (56)) with the above 
function derived from the hyperbolic sine. Letting the subscripts T arid S 
stand for the tangent and sine solutions, we see from equations (54) and 
(86) that for a given Sq, 


2Ayy - Ay s . 


(87) 
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The maximum inverse length scale L as defined by equation (2) 
occurs at t = 0 for the tangent solution, and has the value 

((L _1 . r ) ) T = 2Ay T tanh Ay T . ( 8 ' 

vv t£/max y T J T T 

For the sine solution, the maximum occurs at t = 1 , and has the value 


((L 


-1 


t £.)max) s = Ay s tanh Ay<. . 


(89) 


Thus, for Sq large enough so that tanh Ay^ - tanh Ay<. - 1, the two 
solutions have the same maximum inverse length scale. The minimum slope 
(d£/dt) min. occurs at t = 1 for both solutions. The results are 


( (d?/dt)min) T = 


tanh Ay, 


T Ay, 


and 


tanh Ay<. 

((d5/dt)min) s = ^ 


For large Sq, we thus obtain 


(90) 

(91) 


( (d£/dt)min)<- = l/2( (d£/dt)min) T . 


(92) 


The one-sided function derived from the hyperbolic tangent thus has more 
points at the unclustered end ( t = 1) than the one derived from the 
hyperbolic sine, for identical clustering at t = 0. The difference is 
due to the fact that the zero inflection point occurs at t = 1 for the 
first, and t = 0 for the second. The particular application would 
determine which of these two is preferable. 
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V. CONCLUDING REMARKS 


In this work it has been assumed that the metrics and Jacobians that 
arise in the transformed equations are calculated by finite differences. 
If the equation r(t) of the £ -coordinate curve is known analytically, 
and the transformation t(£) is also given analytically, then the metrics 
and Jacobians can be analytically determined from the derivatives 
dr/dt and dt/d£. The truncation error in the numerical calculation will 
then be due solely to the finite-difference approximations to the 
derivatives of <f> with respect to £ . When <1> varies monotonical ly with t, 
the optimum transformation would be one in which £ varied linearly with 
4> , since this would result in zero truncation errors. 

In order to compare transformations for the numerical and analytic 
treatment of Jacobians and metrics, consider a strictly one-dimensional 
case in which the single unknown u varies monotonically with distance x. 
Assume a highly localized interior region of rapid variation whose 
thickness is proportioned to the small parameterv. A simple example 
of such a solution is 

u ~ tanh (x/v), (93) 

where x = 0 in the localized region. This is actually the steady-state 
solution of Burgers' equation with fixed end conditions. For 

the analytic treatment of Jacobians and metrics, it follows that £(x) 
should be a scaled protion of the hyperbolic tangent. The analysis of 
Appendix E, which assumes a numerical treatment of Jacobians and metrics, 
shows that this choice is completely unsuitable, and instead favors a 
scaled portion of the inverse hyperbolic sine. If the differential 
equation is written so that only derivatives of u appear, then the 
hyperbolic tangent transformation should lead to a numerical solution 
with no truncation errors. But the equation can also be written in a 
form involving derivatives of several functions of u. An example is 
a strong conservation form. In this situation one has several variables 
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(f)(5) to be approximated by finite differences, and no single transformation 
t (5) is optimum for all of them. If vis very small, the hyperbolic tangent 
transformation would put all the interior grid points inside the region 
of rapid variation. There would be no points to resolve the boundary of 
this region. By contrast, the inverse hyperbolic sine transformation 
puts a sufficient number of points outside the region of rapid variation 
to resolve the complete one-dimensional region. 

The above discussion indicates that there are special situations and forms 
of the differential equations for which an analytic treatment of Jacobi ans 
and metrics can provide a desired accuracy with fewer grid points than 
the numerical treatment. These cases appear to be restricted to monotonic 
distributions that can be approximated by simple analytic expressions. For 
general applications of one-dimensional stretching functions, these 
special situations will not be met. It is then best to treat the Jacobians 
and metrics numerically, and use the stretching functions derived in 
Sections III and IV. 

Another assumption in the derivation of the stretching functions is that 
the dimensionless length scale of the localized region of rapid variation 
could be extremely small. This requires the dimensionless slope of the 
transformation dg/dt to be extremely large. If this condition is not 
encountered, and transformation slopes remain of 0(1), then the form of the 
stretching function is not critical. For example, many authors have used 
a scaled exponential as a one-sided stretching function. This is perfectly 
reasonable as long as the one-sided slope s^ is not much larger than one. 

But one can readily show that the maximum inverse length scale is exp(Sg) 
for large Sq. Thus a simple exponential does not yield a suitable 
one-sided stretching function for very large slopes. 

It should be made clear that Roberts was not the first one to use a 
stretching function involving the hyperbolic tangent. Mehta and Lavan 
(ref. 10), in an investigation of flow in a two-dimensional channel with 
a rectangular cavity, used a transformation hased on the hyperbolic 
tangent to transform a semi-infinite region into a finite computational 
region, and to cluster grid points at the corners of the cavity. The 
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maximum dimensionless slope, based on the length of the cavity, had a 
value of .8664 in their calculations. The transformation would have been 
a poor one if they had required a very large slope at the corner, as 
shown in Appendix E. The same authors (ref. 11) used a stretching function 
based on the inverse hyperbolic tangent in a study of the two-dimensional 
flow around an airfoil. A previous transformation had transformed the 
region external to the airfoil into a unit circle. The stretching function 
was necessary to cluster points further near the airfoil surface (to 
capture the boundary layer) as well as the free stream (to overcome the 
stretched grid produced by the first transformation). In their calculation, 
the non-dimensional slopes at the two ends were 5.77 and 27.8. In this 
instance, the use of the hyperbolic tangent was both appropriate and 
necessary, as shown by the analysis of Section III. 

An important criterion in the development of a two-sided stretching 
function is a continuous behavior as s Q and varied from zero to infinity. 
This is necessary to obtain smooth grids constructed algebraically using 
one-dimensional stretching functions. For B > 1 , the required function was 
found to be based on the inverse hyperbolic tangent. At first glance, 
the same function could be used for B < 1 , simply by interchanging £ and t 
in the expression. (This is what was actually done in the earlier stages 
of thiswork). But this would violate the desired continuous behavior 
in the neighborhood of B = 1. The analytic continuation of the inverse 
hyperbolic tangent is the inverse tangent, which differs from the hyperbolic 
tangent. One can actually construct anti -symmetric function which are 
self invertible in this sense, but they do not include the elementary 
functions, and therefore would not be useful as stretching functions. 

Since the inverse tangent is not self invertible, it is necessary to use 
two different representations in calculating the stretching function 
numerically. 

The use of the stretching functions derived in this work requires 
specifying their slopes at one or two points. The values are either 
obtained from matching slopes with another function, or estimating the 
length scale of a localized region where an appropriate dependent variable 
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undergoes rapid variation. Admittedly, in a complex situation, the 
derivation of an appropriate length scale is not easy, and the value to 
assign to the slope can be somewhat arbitrary. Nevertheless, if a 
consistent criterion is used in assigning slope values, useful grids for 
numerical calculations can be generated. A recent example of a complex 
grid which was generated using the general two-sided stretching function 
is found in reference 12. 
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0.2 



FIGURE I TWO-SIDED STRETCHING FUNCTION FOR ^=100 
AND DIFFERENT VALUES OF S, . 
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APPENDIX A 


EVALUATION OF ANTI-SYMMETRIC TWO-SIDED STRETCHING FUNCTIONS 


In this appendix we evaluate several candidates for an anti -symmetric 
two-sided stretching function £(t), where ? and t are normalized variables 
ranging from zero to one. The common slopes at the ends will be 
designated as 


B ■ ft <°> 


dc 

dt 


( 1 ) 


(A. 1 ) 


The truncation error analysis of Section II results in the conditions 

, -1 


ts 


and 


-1 


te 


d 2 t /dt 

d 5 2 / 

= 

d 3 t /dt 

H = 

d 5 3 / 



d 2 c 


dt 

d 

dt : 


( d£,2 

l dt ; 


= 0 ( 1 ) 


<af) 3 - 3L' 2 


t5 


(A. 2) 

=0(1), ( A - 3 ) 


even if B is very large. The corresponding portion of a universal odd 
function w(z) ranges from z Q = -A;/ 2 to z-j = Az/2, where Az is the total 
range, and z is either real or pure imaginary. The stretching function 

is therefore defined as either 

1 

(A. 4) 


w [ Az ( t - j)] 
2 w (az/2) 


° r t 1 - w 

i - 2 i 


2 w (az/2) 


(A. 5) 


Since conditions (A. 2) and (A. 3) are difficult to satisfy only for very 
large B, we restrict the analysis to the case B>1. 
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w = sinz 


If we let z be real, so that Az = Ax, the appropriate function B>1 
is obtained from equation (A. 5) as 

t - i = sin t AX (S' 

2 2 sin (ax/2) • 

Using equation (A.l), we obtain the relation 


„ _ tan (a x/2) 

B " A x/2 

The inverse length scale L"^ takes the form 


L _1 t? = ax | tan [ ax (c - ^)] | , 


and has a maximum value 

(L " 1 t 5 ) max = Ax tan (Ax/2) - 


(A. 6) 


(A. 7) 


(A. 8) 


(A. 9) 


For large B, Ax -* ir, and we obtain 

(A. 10) 

Thus equation (A. 2) is violated, and the function (A. 6) is not suitable. 


< L '\cW * A/2 - 


If we let z be imaginary, so that Az = iAy, the appropriate function for 
B>1 is obtained from equation (A. 4) as 


5 



sinh [Ay(t - Jr)] 
sinh (Ay/2) 


(A. 11) 


Applying equation (A.l), we obtain the relation 


Ay/ 2 

tanh ( a y/2). 


(A. 12) 
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The inverse length scale L ^ is now 

zt > 1 

i -1 = 2 sinhUy/2) | sinh [Ay (t - ~-)] | 

L ££ . _ 

1 + sinh [ Ay(t - ^-)] 

The maximum value of L ^ is readily found to be 

(L ' 1 U ) max = sinh (Ay/2). 

For large B, Ay/2-> B, and we obtain 

(L " 1 U ) max* sinh B ~ e B /2 . (A.15) 

Since the maximum inverse length scale becomes exponentially large, 
the function (A. 11) is completely unsuitable. 


(A. 13) 


(A. 14) 


w = tanz 


The appropriate function for real z is obtained from equation (A. 4) as 


e _ 1 _ tan [ Ax(t - ^-)] 

2 2 tan (Ax/2) . 

The relation of B to Ax is obtained from equation (A.l) as 


(A. 16) 


B = Ax/sin ax. 

The inverse length scale L - ^ is found to be 
1"^ = 2 tan (a x/2) | sin [ax (2t - 1)]. 

For B >tt/ 2, the maximum value of L ^ is given by 

'max ■ 2 ta " 


(A. 17) 


(A. 18) 


(A. 19) 


30 



For large B, Ax ■* it, and equation (A. 17) can be rewritten as 


D 2 sin (Ax/2)cos (A x/2] 
We thus find that 


(7r/2)tan (Ax/2) 


< L ' ■ 


(A. 21) 


and consequently the function (A. 16) is also not suitable. 


The corresponding function for imaginary z is obtained from equation (A. 5) 


* - 


(A. 22) 


The parameter B is now related to Ay by 


sinh Av 


(A. 23) 


The inverse length scale L~ becomes 


L" t £ = 2 Ay | tanh [Ay(£-l/2)] | , 


(A. 24) 


and has the maximum value 


(L ’ t^max = 2Ay tanh (A * /2) • 


(A. 25) 


For large B, equation (A. 23) can be inverted approximately to yield 


Ay -*■ log (2B log B). 


(A. 26) 


(A more accurate relation is derived in Appendix B.) For sufficiently 
large B, although Ay is logarithmically of 0(1), it is large enough for 
tanh (Ay/2)->-l. Consequently we obtain 
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(A. 27) 


('-' Vw 2 '*- 2 109 (2B 109 B> - 


Thus L~ . r remains logarithmically of 0(1), even when B becomes 

tq, 

extremely large. (For example, if B = 1101 , ( L ^)max = 20.) 
Condition (A. 3) requires the determination of L - ^ , which is given 
by 


L - 1 15 = fay 


3 tanl-i [Ay ( ^-l/2)]-l 


1/2 


(A. 28) 


It follows readily for large B that 


( L t^max ~ 


t£^max 


(A. 29) 


and condition (A. 3) is also satisfied. The function (A. 22) is thus a 
suitable candidate for a stretching function. 


w = erfz 


The appropriate function for real z is obtained from equation (A. 5) as 

* - 1' 2 ■ (# - 3o) 


Using equation (A.l), we obtain the relation 


B = /r rerf (Ax/2)e 


Ax 2 /4 


AX 


-1 


The inverse length scale L ^ has the simple form 


= 2Ax 2 1 5- 1/2 1 , 
and has the maximum value 


(A. 31) 


(A. 32) 


< L ' 1 te ) max = ax2 - 


(A. 33) 
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Approximate inversion of equation (A. 31) for large B results in 


Ax 2 + 4 log [2B(log B/tt) 1/2 ] . 


(A. 34) 


We again find that L , f is logarithmically of 0(1), and can also show 
— -1 -1 ^ 

that(L ^Jmax - (L j.g.)max. The error function is not as simple as the 
hyperbolic tangent, and is not invertible. Furthermore, for a given value 
of B, the maximum inverse length scale for the error function is larger 
than that for the hyperbolic tangent. Thus the function (A. 22) is the best 
candidate for a simple, anti -symmetric, two-sided stretching function. 



33 



APPENDIX B 

INVERSION OF y = sinh x/x 


Given the function 
y = sinh x/x. 


(B.l) 


we seek an approximate inverse function x = f(y) for y > 1 . The accuracy 
of the approximation will be measured by calculating the fractional 
error in y, given by 

error (y) . -1. (B.2) 

An expansion for f(y) near y = 1 is readily obtained, but it has a 
small radius of convergence. An asymptotic expansion for large y is 
found to converge very slowly. The range of y for which a high accuracy 
is required extends to y -100. A maximum absolute error of .0005 is 
probably sufficient for the numerical applications of this function. 
Hopefully, these criteria can be satisfied by matching appropriate 
expansions for small and large y with the exact solution at some 
intermediate point y^ . Since the asymptotic expansion has slow 
convergence, it will also be matched to the exact solution at another 
point y 2 - An outline of the derivation of the two expansions will 
now be presented. 


For given values of x-j and x 2 , the corresponding values of y 1 , y 2 , 
(dx/dy) 1 ,(d 2 x/dy 2 )-j ,and (dx/dy) 2 are obtained from equation (B.l) 
by substitution and implicit differentiation. We first consider the 

o 

expansion for y < y-j . Since y-1 + x /6 for small x, f(y)^/6(y-l ) for 
y near one. If we introduce the variable 

y = y - 1 , ( B * 3 ) 

we can write an expansion for x in the form 
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The coefficients A-j and A^ are determined by inverting the expansion of 
equation (B.1) for small x and equating coefficients of y. The results 
are 


A 


1 


= = -.15 


(B.5) 


and 


3x107 

2 32x175 . 


(B.6) 


2 2 

By differentiating equation (B.4), we can set x, dx/dy, and d x/dy equal 
to their known exact values at y^ . The resultant three simultaneous linear 
algebraic equations are readily solved for A^, A^ and A^. 


In order to obtain the expansion for y>y-j,we must first determine the 
leading terms in the asymptotic expansion for large x. Equation (B.l) 
can be rewritten as 

x = sinh" 1 (xy). (b.7) 

Since sinh -1 z ~ log (2z) asymptotically for large z, we obtain the 
asymptotic expression 

x ~ v + log (2x) 

where 

v = log y. 


(B.8) 


(B.9) 
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Substituting the zeroth order approximation 


*<»> 


v 


(B.10) 


into equation (B.8), we obtain the next approximation 


x^ 1 K v + log (2 v) • 

The further substitution of x'^ into equation (B.p>) results in 
x (2 L v + log [2 v( 1 + log(2v)/v)] 

~v + log(2v) + log (2 v)/v 


(B.ll) 


or x* 2 Lv+ (1 + 1/v) log (2v). (B.12) 

The asymptotic approximation (B.12) agrees very well with the exact 
solution for y^l, giving a maximum absolute error of ,02 occurring at 
y~ 200. The next higher order approximation gives poorer results in 
our range of interest. In order to obtain better accuracy, and join 
the solution with the expansion for small y, we add a polynomial in 
inverse powers of y. Since we will also match exact conditions at 
y it is more reasonable to define instead the variable 


w = 1/y - l/y 2 . 

The expansion for x is thus taken to be 

x = v + (1+1/v) log (2v) + B q + BjW + B 2 w 2 + B 3 w 3 + B 4 w 4 . (B. 14) 

The coefficients B q and B-j are determined by equating x and dx/dy to 
their known exact values at y^. Similarly, the coefficients B 2 , B^ and 
B^ are obtained by equating x, dx/dy, and d 2 x/dy 2 to their known values at 

y r 
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For fixed and y^, as y increase above one, the error determined by 
equations (B.2) and(B.4) becomes positive, reaches a maximum, and 
decreases to zero at y^ . Using equation (B.14), we find that the error 
becomes positive again beyond y^ , reaches a maximum, decreases through 
zero, reaches a minimum (which is negative), and increases to zero 
at y^. Beyond y^ the error grows monotonically negative. Thus the 
absolute error reaches a local maximum three times in the range 
1 <y<y 2 . The optimum choice of y-j and y^ is that which makes these 
three maxima equal. By a trial and error procedure, this condition was 
found for the values x-j = 2.722567 and x^ = 6.05012, corresponding to 
y.j = 2.7829681 and y^ = 35.053980. The corresponding maximum absolute 
error is .000267732, which is more than sufficient. This value is again 
reached at y = 69.64. The magnitude of the error increases beyond this, 
reaching .0006 at y~100 and .00083 at y = 120.5. Thus, the desired 
criteria are essentially satisfied by the two approximations. The 
numerical values for the final coefficients are given by equations 
(65) and (67). 
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APPENDIX C 


INVERSION OF y = sin x/x 


Given the function 


y = sin x/x, 


(C.l) 


we seek an approximate inverse function x = g(y) in the range 0<y<l. 
The accuracy of the approximation will be measured by calculating the 
fractional error in y, given by 


error (y) = 



- 1 . 


(C.2) 


The approximation is derived by matching appropriate expansions for y 
near zero and y near one with the exact solution at some intermediate 
point y-j . The derivation is briefly outlined below. 

2 2 

For a given x-j , the corresponding values of y-j , (dx/dy)^ ,and (d x/dy )-j , 
are obtained from equations (C.l) by substitution and implicit 
differentiation. The approximation for y<y-| utilizes the expansion of 
equation (C.l) near y = 0, as a polynomial in powers of tt-x. Based on 
the inversion of this expansion, we take for an approximation the 
expression 


x = tt [ 1 - y+y 2 -(l+ tt 2 /6) y 3 + B 4 y 4 + B 5 y 5 + B 6 y 6 ]. (C.3) 

2 2 

By differentiating equation (C.3), we can set x, dx/dy, and d x/dy equal 
to their known exact values at y^ The resultant three simultaneous 
linear algebraic equations are readily solved for B^, B^ and B g . 


The approximation for y>y 1 similarly utilizes the expansion of equation 
(C.l) near y = 1 . Introducing the variable 


y = 1 - y> 


(C.4) 
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we obtain an approximation in the form 


x =/6y (1 + A-jY + A 2 y 2 + A 3 y 3 + A 4 y 4 + A 5 y 5 ). (C.5) 

The coefficients A-j and A 2 obtained from the expansion near y = 1 are 

A-j = .15 
and 

a 3x107 
rt 2 32x175 . 

The coefficients A 3> A 4 and A 5 are obtained by equating x, dx/dy, and 
d 2 x/dy 2 to their known values at y, . 

The error defined by equation (C.2) reaches a local maximum between y = 0 
and y = y 1 , and a local minimum (which is negative) between y = y-j and 
y = 1. The optimum choice of y-j will make the magnitudes of these extrema 
equal. This was found by trial and error to be produced by x-j = 2.428464 
(corresponding to y-j = .26938972). The maximum absolute error is .000197170, 
which is more than sufficient for numerical calculations. The numerical 
values for the final coefficients are given by equations (71) and (72). 


(C.6) 


(C.7) 
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APPENDIX D 


SOLUTION CHARACTERISTICS OF THE TWO-SIDED STRETCHING FUNCTION 

DERIVED FROM w = tanz 

In this appendix we study the two-sided stretching function derived 
from 


w = tanz, (d.i) 

where z is the complex variable 

z = x + iy- (D.2) 

If Zq and z-j are the two end points of the solution in the z-plane 
which define the range 

Az = z i ‘ z o, (D - 3) 

it is shown in Section III that the appropriate variables for the 
stretching function are defined as 

5 = (z - z Q )/Az (D.4) 


and 


t = (tanz - tanz 0 )/(tanz-| - tanz Q ) > (D.5) 

while the governing parameters A and B are related to Zg and z-j , through 
B = sinAz/Az (D.6) 


and 


A = 


cosz 

cosz 


0 

1 


(D.7) 
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Using equation (D.3), A can also be expressed in the form 


1/A = cosAz - tanzg sinAz. 


(D.8) 


It is also useful to express tanZg as a complex number in the form 


sin 2 x q + isinh 2y Q _ ( D>9 ) 

tanz 0 ~ cos 2 x q + cosh 2y Q 

We will first study the real representations of the solution for 
various ranges of A and B. If B<1, it follows from equation (D.6) that 
Az is real. Using equations (D.8), (D.9), and (D.4), one can prove 
that tanzp, z Q , and z are all real. The stretching function is thus a 
scaled portion of 


w = tanx . 


(D.10) 


This representation breaks down as B + l, since Az + 0. For A<1, it follows 
from equation (D.8) and (D.7) that z+±tt/2. We can obtain the correct 
representation in the limit by introducing 


x = x + tt/2 , 


(D.ll ) 


where x > 0. Since tan (x±n/2)= -cot x--l/x as x + 0, it follows that 
for B = 1 the stretching function is a scaled portion of one branch of 
the rectangular hyperbola 


w = -1/x. 


(D.12) 


If B>1, it follows from equation (D.6) that Az is pure imaginary, 
i.e. Ax = 0. From equation (D.4) it then follows that x = Xq, so that 
the real part of z is constant. Since equation (D.8) required tanzQ 
to be pure imaginary, one can deduce from equation (D.9) that x =-tt/2, 0, 
or +tt/2. If A is sufficiently close to one, the appropriate solution 
is x = 0, and equation (D.l) becomes 


w = itanhy . 


(D.13) 
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Thus the stretching function is a scaled portion of the hyperbolic 
tangent. As A departs from one, this representation must break down. 

One can easily show from equation (D.8) that as A-*exp(-Ay) , yQ^ + 00 
and as A-*exp(Ay), y Q ■+ - « . Since |y|->°°in these two cases, 
equation (D.13) takes the limiting form 

w = + i(l - 2e' 2|y| ). (0.14) 

The stretching function then becomes a scaled portion of an exponential. 
This limit is reached when B attains the critical value 

B* = sinh j log A] 

log A 

For A near one, equation { D . 1 5 ) has the approximate form 
log B* = 1/6 (log A) 2 , 

while for A»1 and 1 / A >> 1 we obtain the asymtotic form 

log B* ->• 1 1 og A| - log ( 2 1 1 og Aj). (0.17) 

For a fixed A<1, as B decreases below B* the appropriate solution for 
z is 


(D.15) 


(D.16) 


z = -J+iy_ (0.18) 

The sign of the correct branch follows from equation (D.ll) by taking 
the limit B-+1. Equation (D.l) then becomes 

w = icothy. (0.19) 

Thus, for 1 <B<B*, the stretching function is a scaled portion of one 
branch of the hyperbolic cotangent. The boundaries for the various 
real representations of the solution are shown in figure D.l, which is 
a log-log plot of the B vs A plane. 
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Another interesting property is to determine the conditions under 
which the solution curve contains an inflection point. This is certainly 
true for the anti -symmetric curve (A = 1). As A departs from one, the 
inflection point moves towards one end, and eventually could disappear. 
The critical point is reached when either 2 y 0 or Zq = 0. From equations 
(D.3) and (D.7), it follows that the former condition occurs when 


A = cosAz . 


(D.20) 


The inflection point thus occurs at the right end of the curve when B 
attains the critical value 


■/ 


A 2 - 1 / cosh' 1 A (A > 1 ) 


(D. 21 a) 


and 



/ cos 


A 


(A< 1). 


(D.21 b) 


For A near one, equation (D.21) has the approximate form 


log B = 1/3 log A, 


(D . 22 ) 


while for A»1 we obtain from equation (D.21a) the asymtotic form 


log B +log 2A - log (2 log 2A). 


(D.23) 


It also follows from equation ( D . 21 b) that A-+0, B attains the 
minimum value 


B . = 2/tt . 

min 


(D.24) 


The condition Zq = 0 is found to occur when 


A = 1/cosAz 


(D.25) 
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Thus the inflection point occurs at the left end of the curve when B 
attains the critical value 


B" =/l - 1/(A 2 ) / cos -1 (1/A) (A > 1 ) 


(D. 26a) 


and 


B" =/l/(A 2 ) - 1 / cosh' 1 (1/A) (A < 1 ). 

For A near one, equation (D.26) has the approximate form 


(D. 26b) 


log B~ = -1/3 log A , 


(D.27) 


while for A«l, we obtain from equation (D.26b) the asymtotic form 


log B' -*■ -log (A/2) - log (-2 log(A/2)). v 

Condition (D.24) is again reached when A-*°°. 

From the above relations, one can conclude that the solution curve 
contains an inflection point unless B lies between B + and B". 

Alternatively, an inflection point will be present if 1/coshAy < A < coshAy 

for B > 1 , and cosAx < A < 1/cosAx for B<1. An inflection point must 

always occur if B < 2 /tt . The regions of the B-A plane where the solution curve 

has inflection points are also shown in figure D.l. Note that these 

regions do not cover completely the regions where solutions are based on 

the tangent and hyperbolic tangent. Thus there are narrow ranges for 

which a solution curve based on the tangent and hyperbolic tangent does 

not contain an inflection point. 

Still another way to study the solution is to follow its path in the 
complex z-plane. This is best done by finding the locus of the solution 
for constant B, as A increases from zero to infinity. It is easy to 
show that as A-*0, x^->-tt/ 2 for all solutions. Similarly, as A->-», 
z^-tt/2 for all solutions. The solution paths are shown in figure D.2, 


44 



with a dashed line for B<1 and solid line for B>1. The details of 
the solution histories are presented below. 

B < 1 


Since z is real, the solution remains on the real axis. When A = 0, 
Xg = -tt/2, and x^ = -ir/2 + Ax. (If B > 2/it , x-j <0, and the solution 
has no inflection point). As A increases, both x Q and x-j increase. 
(If B > 2/u, x-j = 0 and an inflection point appears when A = cosAx) . 
When A = 1, x Q = -Ax/2, x-j = Ax/2, and we have the anti -symmetric 
solution, (if B > 2/n, x Q = 0 and the inflection point disappears 
when A = 1/cosAx). When A reaches +°°, x-| = tt/2, and Xg = it/2 - Ax. 


B > 1 


With Az = iAy, the solution at A = 0 has Zg = -ir/2, and z-j = -tt/2 + iAy. 
As A increases, the solution stays on the line x = -tt/2, and both y Q 
and y-| increase. When A = exp(-Ay), yg = y-| = +°°. As A increases 
beyond exp(-Ay), the solution returns along the imaginary axis, and 
both yg and y^ decrease from +°°. When A = 1/coshAy, y^ = 0, and an 
inflection point first appears. When A = 1 , y^ = -Ay/2, y^ = Ay/2, and 
we have the anti -symmetric solution. When A = coshAy, y^ = 0, and the 
inflection point disappears. When A = exp(Ay), yg = y^ = -°°. As A 
increases beyond exp(Ay), the solution returns along the line x = tt/2, 
and both yg and y^ increase from -». When A reaches +°°, z^ = tt/2, and 
z Q = tt/2 - iAy. 


B = 1 


When A = 0, Zg = z 1 = -tt/2. As A increases, the solution remains at that 
point. When A = 1, z Q and z^ jump to the origin. When A>1, z Q and z-j 
jump to +tt/2, where they remain as A increases to +<». Since B = 1 is 
the boundary between B>1 and B<1, the solution can occur only at these 
points where the B>1 and B<1 solutions intersect. As shown in figure 
D.2, the solid line (B>1) and dashed line (B < 1 ) intersect precisely 
at the three points defined above. 
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FIGURE D.l SOLUTION CHARACTERISTICS OF TWO- 

SIDED STRETCHING FUNCTION DERIVED 
FROM W * TAN Z 



REGION WHERE SOLUTION 
HAS INFLECTION POINT 




FIGURE D.2 : LOCUS OF SOLUTION IN Z- PLANE 


X 
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APPENDIX E 


EVALUATION OF ANTI-SYMMETRIC INTERIOR STRETCHING FUNCTIONS 

In this appendix we evaluate several candidates for an anti -symmetric 
interior stretching function C(t), where 5 and t are normalized 
variables ranging from zero to one. The slope at the midpoint will 
be designated as 

e-af 0/2). (E 

and only the case B>1 will be considered. The criteria to be 
satisfied are that the inverse length scales L~^ and (defined 
in Appendix A, equations (A. 2) and (A. 3)) be of order one, even if B 
is very large. We seek a solution that is a scaled portion of a 
universal odd function w(z). In terms of the range Az, expressions 
for Cor t are given by equations (A. 4) or (A. 5). 

w = tanz 


If we let z be real, so that az = Ax, the appropriate function for B>1 
is obtained from equation (A. 5) as 

t _ 1/2 = tan IaxU-1/2)] . 

1 x/c 2 tan (Ax/2) 

Using equation (B.l) we obtain the relation 

D _ tan (Ax/2) . 

B ' ~Ax/2 

The inverse length scale L”"*^ takes the form 

= 2Ax | tan [Ax(£-l/2 )] | , 
and has a maximum value 

^'Vmax = 2AX tan (Ax/2) - 


(E.3) 

(E.4) 

(E.5) 
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For large B, Ax tt and we obtain 


0 - - ' ) 


t£'max 


■, 2 b. 


(E.6) 


Since the inverse length scale is proportioned to B, the function (E.2) 
is not suitable. 


If we let z be imaginary, so that Az = iAy, the appropriate function 
for B>1 is obtained from equation (A. 4) as 


€ - 1/2 


tanh [Ay (t- 1/2)] 

2 tanh (Ay/2) . 


Applying equation (E.l), we obtain the relation 


(E.7) 


r = Ay/ 2 (E-8) 

tanh (Ay/ 2) . 

The inverse length scale L is now 

L -1 t £= 2 tanh (Ay/2) |sinh [Ay(2t-l)]|, 
and has a maximum value 

(L'^max = 2 tanh ( A y/ 2 ) sinh (E-10) 

For large B, Ay/2-*-B, tanh(/\y/2) +1 , and we obtain 

"-■’tsU* 2 s1nh 2B - e2B - (E - n > 

Since the inverse length scale grows exponentially with B, the function 
(E.7) is completely unsuitable. 


w = sinz 


The appropriate function for real z is obtained from equation (A. 4) 
as 
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r 


1 Z 2 _ sin [A 
' 2 sin 


(E.12) 


t-i/2)3 

Ax/2) 


The relation of B to Ax is obtained from equation (E.l) as 


r = Ax/2 

sin (Ax/2) 


(E.13) 


Since Ax < tt if function (E.12) is to remain monotonic, it follows from 
equation (E.13) that BItt/2. Therefore a monotonic stretching function 
for large B cannot be obtained using function (E.12). 


The corresponding function for imaginary z is obtained from equation (A. 4) 
as 


t - I/? = sinh [Ay (g- 1/2)] 

L w 2 sinh (Ay/2) . 

The parameter B is now related to Ay by 


r = sinh(A.y/2) 
b " Ay/2 

The inverse length scale L - ^ becomes 

L _1 t? = Ay tanh [Ay(£-l/2)], 
and has the maximum value 


(E.14) 


(E.l 5) 


(E.16) 


( L " 1 tc ) max = Ay tanh (Ay/2) - 

For large B, equation (E.l 5) can be inverted approximately to yield 


(E.17) 


Ay -* 2 log (2B log B) . 

For sufficiently large B, Ay is large enough for tanh ( Ay/ 2 ) -»- 1 . 
Consequently we obtain 


U'J^max 1 °9 (2B log B). 


(E.19) 
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Thus L" 1 ., remains logarithmically of 0(1), even for very large B. The 
ts __1 

inverse length scale L ^ becomes 

r ” 1 = Ay . (E.20) 

It follows that for large B 

^ t£^max ~ tf^rnax. (E.21) 


Since both inverse length scales remain of 0(1), the function (E.14) is 
a suitable candidate for a stretching function. 
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